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LINEAR  AND  NONLINEAR  PERSPECTIVES  OF  FORECAST  ERROR 


ESTIMATE  USING  THE  FIRST  PASSAGE  TIME 

Peter  C.  Chu  and  Leonid  M.  Ivanov 
Naval  Postgraduate  School,  Monterey  California 


1.  INTRODUCTION 

Traditionally,  the  prediction  skill  of 
atmospheric  models  is  verified  through  small 
amplitude  stability  analysis.  The  Lyapunov 
exponent  (LE)  and  singular  vector  (SV) 
decomposition  methods  are  the  two  popular 
approaches  (e.g.  Lorenz,  1984,  Dalcher  and 
Kalnay,  1987;  Farrell  and  loannou,  1996; 
Vannitsem  and  Nicolis,  1997  and  others).  The 
model  stability  is  defined  as  sensitivity  to  small 
errors  in  initial  conditions  (the  first  kind  of 
predictability)  and  measured  by  an  e-folding  time 
computed  from  the  leading  LE  (or  SV). 

However,  finite  amplitude  errors  exist  in  many 
practical  cases  for  example,  in  medium-range 
predictions  (Vukicevic,  1991;  Barkmeijer,  1996  ) 
or  “imperfect”  models  (Palmer,  2001  and 
others).  Thus,  the  analysis  should  be  updated. 
Herein,  one  of  possible  approaches  is  the 
probabilistic  analysis  of  forecast  error  dynamics 
(e.g.,  Benzi  and  Carnavale,  1989,  Nicolis  1992; 
Ehrendorfer  1994  ;  Moltineli  and  Corti,  1998  and 
others  ). 

Naturally,  the  knowledge  of  the  probability 
density  function  (PDF)  of  error  allows  to  have 
the  full  statistical  description  of  its  dynamics. 
However,  herein,  even  one  -dimension  error 
dynamics  needs  to  be  studied  by  numerical 
methods  (Ehrendorfer  1994  ;  Nicolis  1992). 

The  probabilistic  approach  can  be 
simplified  if  we  only  determine  the  time  when 
the  model  prediction  skill  is  lost,  i.e.  forecast 
error  became  larger  than  a  given  tolerance 
level  (f)  .  Then,  the  first  passage  time  can  be 
used  as  the  measure  of  prediction  skill  (Ivanov 
et  al.,  1994;  Ivanov  and  Margolina,  1999) 

In  the  present  paper  we  take  such  an 
approach  to  illustrate  its  usefulness 
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for  analyzing  predictability  skill  of  atmospheric 
and  oceanographic  models  and  further 
understanding  how  the  finite  amplitude  error 
affects  low-order  characteristics  of  the  model 
prediction  skill. 

To  do  so  we  use  the  modified  self-consistent 
model  (Nicolis,  1992)  for  error  propagation  and  the 
Princeton  Oceanographic  Model  (POM)  for 
shallow  water  circulation  in  a  semi-closed  basin. 

2.  THE  FIRST  PASSAGE  TIME 

We  a  priori  assume  that  the  model  dynamics 
can  be  described  as  the  stochastic  process  in  an 
N-dimensional  dynamical  system  and  that  the 
characteristics  of  prediction  model  skill  depending 
on  different  physical  factors  are  stochastic.  For  a 
stochastic  process,  the  first  passage  time  is 
defined  as  the  time  when  the  process,  starting 
from  a  given  point  (position),  reaches  a 
predetermined  level  for  the  first  time  (Gardiner, 
1985).  It  was  suggested  by  Ivanov  et  al,  1994, 
Ivanov  and  Margolina,  1999  to  use  the  first 
passage  time  as  one  of  measures  of  the  model 
prediction  skill. 

There  are  several  reasons  to  do  so.  First,  the 
first  passage  time  characterizes  both  linear  and 
nonlinear  perspectives  of  forecast  error.  In  the 
case  of  small  forecast  error  and  for  non 
stochastically  forcing  dynamics  the  mean  passage 
time  coincides  with  the  e-folding  time.  Second,  the 
statistical  moments  of  the  first  passage  time  can 
be  easily  calculated  by  an  iterative  technique 
developed  in  Ivanov  et  al.,  1999.  Third,  in  many 
practical  applications  such  an  approach  allows  to 
get  analytical  estimations  of  the  model  prediction 
skill.  Fourth,  this  approach  is  universal  for  studying 
both  the  first  and  second  kinds  of  predictability. 

Because  the  first  passage  time  is  stochastic 
we  can  introduce  a  special  function  characterizing 
the  model  prediction  skill.  This  is  the  probability 

density  of  prediction  (PDP)  P(t0,^°  ,t  -t0) 
which  is  the  probability  that  the  value  of 
forecasting  error  t,  satisfies  the  inequality 
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<  £  for  the  time  period  (t  —  t0)  dt 

when  it  equals  4  at  t  =  t0.  Herein,  A  is  the 

scaling  matrix,  the  bracket  denotes  the  ensemble 
average  over  realizations  generated  by  stochastic 
forcing  and/or  uncertain  initial  condition.  The 

surface  S  determined  by  =  £ 2  is  the 

ellipsoid  limiting  the  domain  of  possible 
variations  of  forecast  error. 

The  PDP  doesn’t  coincide  with  usual 
probability  density  function  and  satisfies  the 
Pontryagin  -  Kolmogorov  equation  (PKE) 
(Pontryagin  et  al.,  1968,  Ivanov  etal.,  1994): 

LP(t0,$°,t-t0)  =  0  (2.1) 


where  L  =  ^--unV„  - Z)Iffl V (I V ,  V„  =-^- 

df  d£„ 


u,,’Xnmare  determined  from  an  evolution  law  and 

stochastic  forcing  of  hydrodynamic  model  (Ivanov 
etal.,  1994),  n,m=1,...,N. 

PDP  should  satisfy  the  following  initial  and 
boundary  conditions: 

P(t0,Z3°,0)  =  l, 


’t  ^o)  —  0, 

(losing  the  model  skill  when  forecast  error 
reaches  S), 


(impossibility  of  error  to  be  less  than  some  level, 
determined  by  thermal  or  another  kinds  of  noise 
distorted  the  model  dynamics). 

Note  both  boundary  conditions  can  be 
simultaneously  applied  on  different  parts  of  S. 

The  probability  density  of  prediction  can  be 
effectively  calculated  by  numerical  methods,  for 
example,  through  the  ensemble  prediction 
technique  (EPT)  or  by  the  special  iteration 
approach  developed  in  Ivanov  e  t  al.,  1999. 

The  PDP  plays  a  key  role  for  the  estimation 
of  the  model  prediction  skill.  Its  knowledge  allows 
to  determine  the  moment  of  first  passage  time  as 


=  (2.2) 

^0 

A  low-order  characteristics  of  the  prediction  skill 
are  the  mean  and  variance  of  predictability 

time  (r)  =  fj  and  (St2  =  f2  -  T2  respectively. 

Here,  the  brackets  denote  the  ensemble  average 
over  realizations  generated  by  stochastic  forcing. 


The  prediction  skill  for  atmosphere  model  is 
currently  verified  using  the  perfect  model  concept 
with  uncertainty  appearing  only  in  the  initial 
conditions.  In  our  approach,  the  uncertainty  is 
easily  taken  into  account  through  additional 
averaging  of  (2.2)  with  respect  to  an  ensemble  of 
initial  error. 


3.  SELF-CONSISTENT  ERROR  FORECAST 
MODEL 

Nicolis  (1992)  estimated  the  prediction  skill 
of  Lorenz  (1985)  three-component  atmospheric 
model  through  the  self-consistent  error 
prediction  model.  This  model  was  obtained 
through  projection  of  full  error  model  onto  the 
unstable  manifold  and  replacement  of  the 
reference  circulation  to  special  stochastic  term. 
The  model  is  written  by 

where  <z  [0,  °°)  is  non-dimensional  amplitude 
of  error,  <7  and  g,  are  generally  time- 
independent  parameters  whose  properties 
depend  on  the  underlying  model  attractor, 

v(t)  =0,  v(t)v(t')  =q2S{t-t'). 

We  modified  the  self-consistent  model 
because  in  the  real  world  any  dynamics  is 
distorted  by  at  least  the  thermal  (molecular) 

noise  and  value  of  initial  error  4  cannot  be 
less  than  the  amplitude  of  the  thermal  noise  4 . 

Let  us  account  this  noise  and  introduce  the 
special  scale  4  =  ((i ,herm))' 1 2 ,  Itherm  is  the 
intensity  of  the  thermal  noise.  The  contribution 
of  4  to  the  predictability  time  can  very  small  if 

4  » 4  ■  However,  use  of  4  leads  to  a 
correct  asymptotic  of  the  predictability  time  for 
4  0 .  The  introduced  scale  also  can  be 

interpreted  as  the  characteristic  scale  of 
unresolved  sub-grid  motions. 

Although  the  model  (3.1)  is  simple  it 
demonstrates  several  dynamical  regimes  of 
error  behavior.  The  non-dimensional  parameter 

/l  =  2 G t  q2  is  a  threshold  between  different 
dynamical  regimes  of  error  evolution. 

For  A  >  1 , 4  «  1  ar|d  g  «  1  we 


obtained  the  e-folding  time  written  as 

/  \  2  ,  ,  £ . 

<3-2> 

q  (A  1)  4 


for  the  critical  point  A  =  1 
the  dependence  for  the  predictability  time 
differs  from  the  e  -folding  time 


M  =  4  taC4)k«f ")  -  2  ln(4 )]  (3.3) 

r  # 

which  demonstrates  that  reduction  of  any  noise 
intensity  increases  the  predictability  time  but 
only  on  the  order  of  ~  In (Sth )  . 


Lower  the  critical  point  (A>  1)  the 
predictability  time  can  be  determined  as 


T)  =  - 


1 


(l-A)q2 


mp-np 


(3.4) 


where  B( — )  = 


T 

1 

(f° 

y-/i 

x(. 

-A 

.0= 


q+i) 

&i-/i 

Sth 


The  effect  of  nonlinear  terms  in  (3.1)  appears 
in  the  power  terms  in  (3.2)-(3.4),  which  reduce 
the  predictability  time  and  increasing  it  variance. 
Detailed  results  are  presented  in  Chu  et  al. , 
2001a. 


4.  PREDICTION  SKILL  OF  POM 


For  chosen  mean  wind  stress 
Tw  =  10~2  m2s~2  and  basin  geometry  (a 
rectangular  basin  with  the  geometrical  sizes  of 
1 000 km  x  1 050 km  x  2 km  ,  plane  bottom  and  a 
single  open  boundary)  we  found  that  for  50-60 
days  the  flow  starting  from  the  state  of  rest 
develops  in  the  quasi-stationary  inferred 
circulation.  We  ignored  the  horizontal  diffusion, 
and  only  the  bottom  friction  was  considered. 
The  details  of  numerical  experiments  can  be 
found  in  Chu  etal. ,2001  b. 

The  non-stationary  circulation  developed 
within  a  two-month  period  was  chosen  as  the 
reference  solution  for  the  study  of  the  POM 
stability.  We  focused  only  on  the  short  and 
intermediate  predictions  limited  by  2-2.5 
months.  The  reference  solution  were  distorted 
by  stochastic  errors  of  initial  conditions, 
stochastic  wind  and  normal  velocity  along  the 
open  boundary. 

We  calculated  numerically  by  EPT  and 
analytically  from  PKE  the  probability  density  of 
prediction,  the  predictability  time  and  its 
variance. 

The  mean  and  the  most  probable  times,  the 
variance  of  the  mean  time  are  the  local 
characteristics  of  the  probabilistic  model  stability 


in  the  perturbation  geostrophic  stream  function 
L,  norm  . 

In  order  to  determine  the  prediction  skill, 
the  POM  was  reduced  to  a  simpler  stochastic 
dynamical  system.  The  system  is  a  set  of 
stochastic  ordinary  differential  equations  whose 
parameters  are  embedded  in  probabilistic 
properties  of  the  full  model. 

We  have  developed  the  set  of  orthogonal 
basis  functions  (normal  modes)  for  semi-closed 
seas.  That  allowed  constructing  a  dynamical 
system  for  the  full  POM. 

To  do  so  we  represent  velocity  vector 
through  the  two  scalar  potentials  W  and  <P 
(Chu  et  al.,  2001  b)  as 

u  =  Vx(k*F)  +  V^  (4.1) 

and  chosen  the  following  open  boundary 
condition  for  the  velocity  potential 

<P  =  0  (4.2) 

Such  an  approach  allows  to  identify  W  as 
the  geostrophic  steam  functions  and  exclude  the 
velocity  potential  from  the  analysis  because  ratio 
of  mean  kinetic  energy  correspond  to  potential 
velocity  and  mean  kinetic  energy  of  irrotational 
motions  was  less  than  0.01. 

Then,  we  introduced  a  new  variables  by 

v'  =  v-\umda  (4.3) 

0 

where  unorm  is  the  normal  velocity  along  the 
open  boundary. 

The  transformation  (4.3)  allowed  to 
reformulate  the  POM  stability  problem  into 
uncertainty  of  normal  velocity  along  an  open 
boundary  as  the  probabilistic  stability  of  a 
dynamical  system  with  multiplicative  stochastic 
forcing  and  to  introduce  the  spectral 
representations 

M 

(4.4) 

m= 1 

where  ^„,is  the  eigenfunctions  of  Laplacian 
operator. 

We  have  found  that  the  10-  dimensional 
dynamical  system  for  Am  approximated  the 

reference  solution  with  the  accuracy  around  1-2 
%.  In  comparing  with  ten  normal  modes  three 
empirical  orthogonal  functions  contents  more 
than  0.1%  of  total  kinetic  energy  of  the 
reference  circulation. 

Then,  PDP  was  calculated  for  the  full 
model,  for  dynamical  systems  of  10th(M=10)and 
100th  (M=100)  orders  by  EPT  and  only  for  a  10- 


dimensional  dynamical  system  from  PKE.  We 
vary  ensemble  size  between  50  and  10000 
members. 

We  have  found  that  the  forecast  error  can 
be  strongly  reduced  on  the  initial  growing  stage. 
The  zero  growth  (or  even  negative)  of  forecast 
error  within  30-day  period  is  obtained  if 
Gaussian  stochastic  noise  with  the  correlation 
radius  equaled  to  300  km  was  added  into  initial 
conditions  and  uncertainty  of  open  boundary 
conditions  was  small.  Reducing  the  correlation 
radius,  i.e.  transiting  to  white  noise,  we 
reanimated  effective  growth  of  initial 
perturbations. 

The  prediction  skill  of  the  POM  is  the  most 
sensitive  to  the  systematic  errors  of  normal 
velocity  along  the  open  boundary.  Stochastic 
perturbations  of  the  normal  velocity  can  often 
result  only  in  super-slow  growth  of  the 
forecasting  error. 

Another  interest  physical  result  is  that  the 
spectrum  of  forecast  error  collapse  fast  into 
low-frequency  domain  of  wave  numbers.  The 
physical  reason  of  this  phenomena  is 
geostrophic  adjuction  of  the  reference 
circulation. 

Three  typical  scenarios  of  the  predictability 
process  have  been  found. 

Scenario  1.  After  a  rapid  growth  of  initial 
perturbations  a  non  Gaussian  probability  density 
function  (PDF)  for  the  first  passage  time  had  a 
symmetric  form  with  narrow  variance  (  Figure  1 

a,d,e;  £  =  10-1).  Such  a  scenario  is  the  most 
probable  for  large  forecast  error.  Herein,  the 
ratio  of  the  mean  predictability  time  to  the  most 
probable  (a)  is  around  1.  The  mean 
predictability  time  calculated  by  EPT,  Singular 
vectors  and  from  PKE  was  equal  to  27.6,  26.9 
and  25.6  days,  respectively.  The  variance  of 
this  time  written  in  the  same  notations  isl  1 1 .7,  - 

and  10.5  day2 .  The  prediction  skill  of  the  model 
could  be  considerably  improved  through  the 
filtration  of  the  fastest  growing  initial 
perturbations  . 

Scenario  2.  We  illustrate  it  through  the 
analysis  of  flow  forced  by  stochastic  wind  with 
the  variance  equaled  to  21  m 2  /  s2  (  Fig. 2  b,e,h  ; 

£  =  10-2 ).  Herein,  the  initial  condition  was 
weakly  perturbed. 

PDF  had  a  non-  symmetric  form  with  a 
long  tail  in  domain  of  the  long-term  prediction, 
df  —  1.3.  The  calculated  mean  and  variance  of 

the  first-passage  time,  (r ^((St2^),  are  1.32 


(0,31);  and  1.31  (0.33),  respectively.  Although 
the  mean  predictability  time  is  not  large,  the 
predictability  time  for  single  realizations  can  be 
very  large. 

We  called  it  “extreme  prediction”,  which 
strongly  depends  on  properties  of  the  model 
attractor  and  weakly  on  uncertainty  of  initial 
conditions. 

Notice  that  such  a  scenario  is  often  take  place 
for  small  values  of  tolerance  £  ,  and  indicates 
a  linear  evolution  of  forecast  error. 
Improvement  of  the  model  prediction  skill 
requires  some  different  approaches  from 
reduction  of  initial  errors. 

Scenario  3.  This  scenario  is  for  the  case 
that  only  uncertainty  of  open  boundary  condition 
exists.  Instead  of  uniqueness,  PDP  had  several 

maximuma,  a  ~  5.2  (Fig.  1  c,  f,  i;  £  =  5.10-2 ). 
Use  of  the  notations  in  scenario  2, 

was  equal  to  10.6  (84.3)  and  12.1 

(88.3),  respectively.  Herein,  we  have  both 
short-term  and  long-term  (extreme)  predictions. 

Notice  that  the  uncertainty  of  initial  conditions 
plays  a  second  role  again.  The  tail  of  PDF  can 
be  directed  along  both  short  and  long-term 
predictions. 

The  principle  property  obtained  from  the 
scenario  is  abnormal  large  value  of  the  variance 
of  the  predictability  time.  Numerical  calculations 
demonstrate  that  the  scenario  3  can  be  realized 
for  both  the  small  and  large  values  of  forecast 
errors. 

5.  THE  EXTREME  PREDICTIONS  IN 
OCEAN 

In  order  to  check  existence  of  extreme 
predictions  found  numerically  in  the  real  world 
we  estimated  the  prediction  skill  one  of  high- 
resolution  circulation  models  of  the  Gulf  of 
Mexico.  This  model  was  developed  by  Kantha 
and  Clayson,  2000. 

The  model  was  integrated  on  15X15  km  grid 
and  forced  by  monthly  wind  and  assimilated 
climatic  boundary  conditions  in  Yakatan  Channel 
and  Strait  of  Florida.  The  model  prediction 
skill  was  estimated  using  difference  between 
trajectories  of  buoys  and  of  synthetic  particles 
deployed  in  numerically  simulated  circulation. 

We  introduced  the  mean  predictability  time 
averaged  over  buoy  trajectories  and  by  all 
buoys.  Number  of  terms  of  such  an  ensemble 
achieved  1600. 


Our  calculations  clearly  show  the  extreme 
predictions.  For  example,  for  50-  km  tolerance 
level,  and  the  most  probable  time  equals  3 
days,  the  extreme  prediction  time  can  reach 
12-15  days  (Fig .2  a,  b). 

The  mean  predictability  time  linearity  grows 
with  increasing  the  value  of  tolerance.  For 
values  of  tolerance  are  more  than  50  km  we 
had  30-50  day  successful  predictions  of  buoy 
trajectories. 
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